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DECUS Program Library Write-up 


DECUS NO. 8-567 


Sly 

SUB 8K-# Jsau 9, 1972 

(1) NAME : EXPO 

Purpose : To evaluate the approximate state 

transition matrix/ its augmented 
transition matrix, 

AT 

e 

where A is a square matrix, and T 
is a sampling period, which appear 
usually in the state variable approach 
to engineering problems* 

X(nT) * A X(0) 

where X(t} is an m x 1 state 
variable vector at state t, X(G) 
is the initial values of X(t), and 
A is in general time dependent state 
transition matrix of dimension m x m v 
if the input forcing function is not 
.identically zero, then A is treated 
as an augmented transition matrix 
which may or may not be singular* as 
long as the approach of solution avoid 
using inverse oS. A fr which is what 
done in EXPO 0 

(2) CALLING SEQUENCE : 

CALL EXPO(M,T,NPTS,DIGIT) 

where; 

Input Data : 

A(X,J),I,J« 1,2, o.o.M ( M at most 7) 

M s ? of row/ # of column/ # of state 
variables used 
T » sampling period 
NPTS * # of calculated samples 
DIGIT ss desired # of significant figures accuracy 

Return Data : 

AX(I,J), I,J*s1,2, .«►*, M the approximate transi¬ 
tion matrix 

(3) ERROR RETURN : 

n FPNT" ERROR AT LOC 03350 will be printed if 












coi?!putar overflow/underflow# 


(4) SPECIAL gONSIP EBAgCOM : 

If machine overflow/underflow» user is urged to re¬ 
submit the program for process with & reduced value 
of T# 

(5) SUBPROGRAM CALLED ; 

-none- 

(6) STORAGE REQUIR ED: 

pages(octal) 

(7) ALGORITHM & R EFERENCES : 

Algorithm i 
AT 

e* A a M + R K 


where i 


atepl 

otepp 


M 


ii T x 

i»0 

i J 

R 

C^£5 

1 % 

i-rT^ 


lisrr 

i ! 

Choose K arbi 

■••reri! 

the 

elements 

{52.. . 

id 


, A u a I s idnentity 
matrix 


remainder matrix 


Determined e by 
e . 9A| 1 -I 

K+2 


where;: the Linear normals defined as 

Tj M > 


A | a max ( L 

1 i 


step3 ’ Find the upper bound of each element in 
the remainder matrix, 

r | « , i 

U] * (K+1) ! “T~ 

step/f ; Compare each w . from stop 1 with the 
uppoer bound of |r. J from step 3, if 
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where d a accuracy at 


least ,} d" eignif'ic&tnt figures is not 
satisfled, then increase K and repeat 
step 1 thru step 4, otherwise the pro* 
cess is finished. 
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(8) LISTING : 

(see attach) 


(9) SAMPLE : 

(see attach) 
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c **** CALCULATION OF APPROXIMATE ” STATE. TRANSITION MATRIX 
C **** JAN. 9, 1972 

SUBROUTINE EXP 0 C M, T , NP TS , DIGIT > 

COMMON A# AX 

DIMENSION AC7,7),AX<7,7),C(7,7>,E<7,7> 

KK»4 

22 DO 10 1=1 >M 

DO 10 J=1,M 
CfI,J)s0. 

AX<I>J>=0. 

10 CONTINUE 

DO 20 I=1,M 
AXCI,I)=1• 

20 C(I,I)=1. 

A(C=1 . 

DO 50 II=1,KK 
DO 30 I =1 , M 
DO 30 J=1 ,M 
S=0. 

DO 25 K = 1 ,M 

25 S=S+CCI,K)*ACK,J) 

ECI,J>=S 
30 CONTINUE 

C ..... K FACTORIALS 

AK»FLOAT<m*AK 
DO 50 I =1 ,M 
DO 50 J=1,M 
C<I,J)=T*E<I,J)/AK 
AXII,J)=AX<I,J)+CCI,J) 

50 CONTINUE 

AMAX=0.2E-38 
DO 70 1=1,M 
S=0. 

DO 65 J=1,M 

65 S=S+ABS(A(I, J) ) 

IFCAMAX-S) 55,70,70 
55 AMAX=S 

70 CONTINUE 

C *** RATIO OF A TERM BY ITS PRECEEDING TERM 
EX=FLOAT(KK+1) 

EPS=AMAX*T/CEX+1. ) 

C *** UPPER BOUND OF EACH ELEMENTS IN REMAINDER MATRIX 
R=(EPS*CEX+l.))**EX/C <AK*EX)*C1.-EPS)) 

DO 80 1=1,M 
DO 80 J=1,M 

C *** ACCURACY UP TO " DIGIT ” SIGNIFICANT FIGURES 
IF(R-ABSCAXCI,J>)*10.**<-DIGIT)) 80,80,85 
85 KK=3* KK/2 

GO TO 22 
80 CONTINUE 

RETURN 
END 
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i<°‘)= <Lr A ^ 
x(p*)= a^~ - k* <“) 







c **** SAMPLE # 3* 

C **** INVERSE LAPLACE TRANSFORM OF A RATIONAL FUNCTION 
C 

C **** A IS AN M BY K STATE TRANSITION MATRIX 
C M = 0 OF STATE VARIABLES 

C NT * TOTAL # OF PTS» TO BE PUNCH OUT 

C TMX ■ MAX. TIME 

C TSC a TIME SAMPLING PERIOD CALC ( IN ••EXPO” ) MUST 

C BE LESS THAN OR EQUAL TO <TMX/<NT-1 )) 

C 

COMMON Ai C 

DIMENSION A(7*7),CC7*7)*X<7) 

HEAD(2*5> M* NT#TMX*TSC*DIGIT 
5 FORMAT(I10/I10/E12.6/E12.6/E12.6) 

NP«TMX/FLOAT<NT-1>/TSC 
NC=NP*NT 
DO 30 1 = 1 *M 
DO 20 J»1#M 

20 HEADC2#10) A(I*J)* 

FINI 

30 CONTINUE 

10 FORMAT(El 0.4) 

C *** READ STATE VARIABLES INITIAL VALUES 
DO 25 1 = 1 *K 

25 R£AD<2*10) X(I) 

C 

CALL EXPO(M*TSC*NC#DIGIT) 

T2=3 • 

DO 320 N=i>NC 

C *** PRINT TIME* TRANSIENT SOLUTION* FIRST DERIVATIVE 
IF(IREM((N-1)/NP)) 316*312*316 
312 .;HITE(2*315) T2#XU)*X(2) 

315 FORMAT(3E13.4) 

316 DO 310 I=1*M 
S=0 . 

DO 303 J=1*M 
300 S=S+CCI*J)*X(J> 

xcn«s 

310 CONTINUE 

T2=T2+TSC 
320 CONTINUE 

STOP 
END 




C **** A IS AN 0" By ft STATE TRANSITION '/A l'HIX 

C 

COMMON A , C 

DI MENS I ON AC / • () , G < ,V/)^(7) 

READC2,5) M, MRT5, I » DIGIT 
5 FORMA I ( 11.0 / T 1 0 /K l > . > 'l l P » (S ) 

DO 30 1*1,0 
DO A-I J-1. ,m 

20 HKADC2,1 A ) AC 1,0) , 

FIM 

30 CONTINUE 

10 F ORMAT < El 0•A) 

C *** READ STATE VARIABLES INITIAL VALUES 
DO 25 1=1 »M 

25 R EA D C 2,10) XCI) 

C 

CALL EXPOC M,T,NPIS,DIGIT) 

T2=0. 

DO 380 N*1,NPTS 

C *** PRINT TIME, TRANSIENT SOLUTION, FIRST DERIVATIVE 

WR ITE C 2,31 5 ) T2,XC1),XC2) 

315 FORMAT C3E13»4> 

DO 310 1=1,M 
S=0 . 

DO 300 J=1,M 
300 S=S+CCI,J)*X(J) 

XCI)=S 

310 CONTINUE 

T2=T2+T 

320 CONTINUE 

STOP 
END 
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c **** SAMPLE * 5 < JILL OUTPUTS TIKE# INPUT# RESPONSE ) 

C CALLS EXPO & INPUT SUBROUTINES 

C **** SOLVINGS 
C X* ■ A*X ♦ n*u 

C WHERE X IS A STATE VARIABLES 
C A IS TRANSITION MATRIX • 

C B ASSOCIATED MATRIX TO MwPUT FORCING FUNCTION- U 

c «•••• 

C MUST CONVERT TO* 

C 

C X* =» A*X WHERE A IS NOW AN AUGMENTED MATRIX 

C 

c ' ■ * 

C **** A IS AN M BY M "AUGMENTED TRANSITION MATRIX" 

C M a * OF STATE VARIABLES 

C NT = TOTAL 0 OF PTS. TO BE PUNCH OUT CPER VARIABLE) 

C TMX = MAX. TIME 

C TSC « TIME SAMPLING PERIOD CALC < IN "EXPO" ) MUST 

C BE LESS THAN OR EQUAL TO <TMX/(KT-l)> 

C 

COMMON A# C 

DIMENSION A(7#7>#C<7#7>#XC?) 

READC2# 5) M# NT#TMX.TSC.DIGIT 
5 FORMAT Ul0/I10/El 2.6/El2.6/El 2.6) 

NP«TMX/FLOAT(NT-I)/TSC 
NC=NP*NT 

DO 30 I «1 #M ? 

DO 20 J»l#M ; • 

20 READ(2#10) AC I# J># 

FINI 

30 CONTINUE 

10 FORMAT < El 0.4) 

C **** READ STATE VARIABLES INITIAL VALUES 
DO 25 I aJ #M 

25 READC2# 10) XU) 

C 

CALL EXPO(M#TSC#NC#DIGIT) 

T8«0. 

DO 320 N«1#NC 

C **** PRINT TIME# INPUT FCN# TIME RESPONSE 
IFUREMC (N~l )/NP) ) 316.312.316 
312 WRITE<2#315) T2#X(M)#XC1) 

31 5 FORMAT < 3E13 • 4 ) 

316 DO 310 1 = 1 #M 

Sa0 . 

DO 300 «Ja l # M 
300 SaS+C(I#J)*XC«J> 

XU >»S 

310 CONTINUE 

T2*T2+TSC 
C 

CALL INPUTCX4.T2) 

C **** SET INPUT AS LAST VARIABLE 
X<M)=X4 

320 CONTINUE 

STOP 
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END 
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INPUT # 3 

ECT) * 103*T# FOH 0<» T <=*l . 

« 0 FOH T> 1• 


SUBROUTINE INPUTCX4>T2) 
IFCT2-1.) 5/5*13 
X4»103.*T3 

RETURN 

X4»0. 

RETURN 

END 
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